Faraday waves in Bose-Einstein condensates 
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Motivated by recent experiments on Faraday waves in Bose-Einstein condensates we investigate 
both analytically and numerically the dynamics of cigar-shaped Bose- condensed gases subject to 
periodic modulation of the strength of the transverse confinement. We offer a fully analytical 
explanation of the observed parametric resonance, based on a Mat hieu- type analysis of the non- 
polynomial Schrodinger equation. The theoretical prediction for the pattern periodicity versus 
the driving frequency is directly compared with the experimental data, yielding good qualitative 
and quantitative agreement between the two. These results are corroborated by direct numerical 
simulations of both the one-dimensional non-polynomial Schrodinger equation and of the fully three- 
dimensional Gross-Pitaevskii equation. 



I. INTRODUCTION 

Pattern formation in driven systems is an important 
direction of current research that influences many fields 
ranging from hydrodynamics to biophysics and from non- 
linear optics to reaction kinetics; see [1] for a comprehen- 
sive review of the topic. 

Some of the oldest and most well-known forms of such 
phenomena are the so-called Faraday patterns, stemming 
from the classical studies of Faraday in 1831 [2], who 
studied the behavior of "groups of particles [placed] upon 
vibrating elastic surfaces" and (in the appendix of his 
much-celebrated paper) the dynamics of "fluids in con- 
tact with vibrating surfaces." Faraday's experiment be- 
came a classical example of pattern formation, whereby 
the uniform state loses its stability against spatially mod- 
ulated waveforms, whose dominant length-scale is deter- 
mined by the intrinsic properties of the system (such as 
the dominant wavelength of the instability) and is only 
weakly dependent on boundary or initial conditions. 

In the past few years, there has been an increasing lit- 
erature about observing phenomena of the above type in 
driven superfluids. Bose-Einstein condensates 0, 0] offer 
perhaps the ideal playground for such experiments, since 
their experimental tunability permits to create such para- 
metric resonance phenomena in a multiplicity of ways. 
One such way is by driving the magnetic trap confining 
the system, as was proposed in Ref. [5|. In the same 
spirit of modulating the confinement of the system to 
observe parametric resonances, the works of Refs. 0, 
motivated by the experiments of Ref. Q , considered a pe- 
riodic modulation of an optical lattice potential; on the 
other hand, the later work of Ref. @ focused on modulat- 
ing in time the confinement potential in a toroidal trap. 
Another recent suggestion was to produce a parametric 
drive by means of periodically modulating the scatter- 
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ing length [10[, which is directly related to the prefactor 
of the effective nonlinearity (due to the mean-field inter- 
particle interactions) of the system. This can be achieved 
by means of Feshbach resonances [ll| . 

These theoretical propositions motivated the very re- 
cent experimental implementation of the Faraday waves 
in Bose-Einstein condensates in the work of Ref. [l2| ]. 
The actual realization of the spatially modulated pat- 
terns arose in a somewhat different way than was pro- 
posed in the above studies, a way which is very close, 
however, to the spirit of the theoretical suggestion of 
Ref. [l3[ (see also the more recent consideration of a sim- 
ilar problem from a quantum point of view in Ref. [lil]). 
In particular, in Ref. [12[, an elongated cigar-shaped 
condensate was used where the transverse, strong con- 
finement directions were periodically modulated in time, 
while the weaker longitudinal direction confinement was 
time-independent. The parametric excitation at the driv- 
ing frequency was recognized as being responsible for ex- 
citing oscillations at half the driving frequency, which 
is the main resonance also observed in Faraday's ex- 
periments. Subsequently, based on this insight and the 
dispersion relation of longitudinal collective modes pre- 
sented in Refs. [H, [l6[, a relation was derived (and con- 
vincingly compared to the experimental results in a quan- 
titative manner) for the resulting pattern periodicity ver- 
sus the transverse driving frequency. 

The aim of the present paper is to provide a complete 
analysis of the instability from first principles and to ob- 
tain a fully analytical prediction that can be used for a 
detailed comparison with the experimental results and 
numerical results obtained from the ID model reduction 
as well as full 3D simulations at the mean-field level. 

The principal feature which allows us to provide a 
detailed quantitative analysis of the system is the fact 
that for cigar-shaped condensates (such as the ones used 
in the experiment of Ref. [12]) there is a quantitatively 
accurate description in the form of the non-polynomial 
Schrodinger equation (NPSE) derived in Ref. [17| (pro- 
vided that the transverse direction stays close to its 
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ground state, which is approximately the case in the ex- 
periments of Ref. |l2|). The remarkable feature of the 
NPSE is that in the resulting partial differential equation 
(PDE) for the longitudinal description of the condensate, 
the transverse frequency enters explicitly and hence pro- 
vides, in this setting, the explicit parametric drive that 
will lead to the observed spatially modulated patterns. 
To elucidate this feature, we will perform a linear sta- 
bility analysis of the spatially uniform states in this ef- 
fectively one-dimensional setting. This naturally leads 
to a Mathieu equation. Then, one can use the theory of 
the Mathieu equation to identify the most unstable mode 
and the wavenumber (and associated wavelength) of its 
spatial periodicity. This, in turn, allows us to obtain a 
fully analytical expression for the pattern periodicity as 
a function of the driving frequency that can be directly 
compared with the experiment. We believe that this ap- 
proach yields both qualitative and quantitative insight 
on the experiment and on the nature of the relevant phe- 
nomena. 

Our presentation will be structured as follows. In Sec- 
tion II, we develop our analytical approach and directly 
compare it with the experimental findings. In Section III, 
we complement our theoretical findings with numerical 
results simulating both the full 3D experimental setting 
and comparing it with results of its ID analog, namely 
the NPSE equation. Finally, in Section IV, we summa- 
rize our findings and present some interesting directions 
for future work. 



II. ANALYTICAL CONSIDERATIONS 

Our starting point will be the standard mean-field 
model for Bose-Einstein condensates (BEC), namely the 
three-dimensional Gross-Pitaevskii (GP) equation 0, 0] 
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where ip denotes the BEC wavefunction (normalized to 
unity), N represents the number of atoms and g = 
47rh 2 a s /m is proportional to the scattering length a s of 
the interatomic interactions; the external potential con- 
fining the atoms is given by: 
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where r 2 



y 2 , uo r is the transverse trapping fre- 



quency, and U(z), in the case of Ref. [12[, represents 
a far weaker parabolic, longitudinal, potential (with a 
magnetic trap frequency less than 5% of the transverse 
frequency) , which will therefore be neglected for the pur- 
poses of the (main analytical portion of the) present 
work. 

Following the approach of Ref. [13], the three- 
dimensional wave function can be decomposed into a ra- 
dial and a longitudinal part as 



where the radial component is taken as 

/ ex P \-r 2 /(2(J 2 (z J t))} 

0(r, t- v(z, t)) = PL J\ y . )} \ (3) 

with a spatially and temporally variable width charac- 
terized by a(z, t). It should be noted here that the radial 
profiles of high- density cigar-shaped condensates (as the 
one in Ref. [12[) are closer to the Thomas-Fermi regime 
than to the Gaussian one. While this introduces an ele- 
ment of approximation to the calculation below, the ad- 
justable nature of the parameter cr(z,t) (which is always 
larger than the transverse oscillator length by a spatially 
dependent factor that depends on the longitudinal wave- 
function; see Ref. fl7| ) and the comparison that we will 
report below between our theory and the physical ex- 
periment render this approximation a reasonable one for 
the purposes of predicting the wavelength of the result- 
ing pattern. We note in passing that it appears to be an 
interesting openproblem to perform a derivation similar 
to that of Ref. [17[ , under the assumption of a transverse 
Thomas- Fermi wavefunction profile. 

Employing the standard variational recipe of Ref. [ItJ , 
and neglecting the longitudinal potential U(z), one ob- 
tains the following effective PDE, the so-called non- 
polynomial Schrodinger equation (NPSE), describing the 
longitudinal wavefunction: 
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In accordance with the experimental setup of Ref. [l2j], 
the transverse frequency is modulated by 



uo r (t) = c<; r ,o • (1 + esin(o;t)), 



(5) 



where uo r ^ is the reference tranverse frequency, and e and 
ou are, respectively, the amplitude and frequency of the 
modulation. Then, the spatially homogeneous solution is 
given by 



fo(t) = Aexp 
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(7) 



t/>(r,t)=<l>(r,ti*(z,t))f(z,t), 



(2) 



and A is a positive constant. The numerical value of A is 
computed from the normalization JJJ \ip(r, z, t)\ 2 dr = 1. 
Computing the integral one has that A = ^/l/2L, where 
we assume that the condensate extends between — L and 
L; i.e., we are assuming here that the condensate is in 
a box rather than in a very weak magnetic trap in the 
longitudinal direction. The validity of this assumption 
for the present phenomenology is verified both a priori 
(due to the very weak nature of the longitudinal con- 
finement in comparison to the much stronger transverse 
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confinement and its modulation) and a posteriori (from 
the comparison with the experimental results). 

Faraday patterns appear in this context due to a mod- 
ulational instability along the (longitudinal) z-axis. To 
examine the modulational stability of uniform patterns 
in the z-direction, we use the ansatz 

f(t) = f (t) [1 + (u(t) + iv(t)) cos(kz)] . (8) 

Inserting this ansatz into Eq. (jlj) and linearizing the en- 



suing equations yields a Mathieu-type equation for the 
perturbation: 



d 2 u 
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o, 
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where a and b are given by Eqs. (jTUJ) and ([TT]) and r 
ujt/2. 
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As is commonly known, Mathieu equation exhibits an 
intricate stability chart comprising tongues of both stable 
and unstable solutions [l8[ • Due to the periodic potential 
a generic solution takes the form u(t) = e ZfJjt g(t), where 
g(t) has the same periodicity as sin(2r) (according to the 
Floquet-Bloch theorem), and \i is a complex exponent 
taken as /i = /ii + i/i2, where both \i\ and /i2 are real 
numbers. Determining the most unstable mode (which 
is the one that is expected to be observed experimentally) 
amounts to finding the oj(k) curve corresponding to the 
critical exponent with the most negative imaginary part. 
While this is usually a complicated task, it can be shown 
that for small positive values of b it amounts to 



a(fc, uj) w 1. 



(12) 



This conclusive property can be seen both numerically 
and analytically. 

Investigating \±2 numerically pjl one finds that (for 
small positive values of b) it consists of a set of symmetri- 
cal lobes centered around a = n 2 , where n is a positive in- 
teger, the one centered around a = 1 being substantially 
larger than the other ones (see Fig.[T]). These lobes corre- 
spond to the unstable regions. Inspecting these lobes, it 
is transparent that the most unstable mode corresponds 
to a « 1. 

One can also argue for the validity of Eq. ([T2|) by an- 
alytical means. There is a class of partly-forgotten ap- 
proximate formulas for \i as a function of a and b stem- 
ming from celestial mechanics (see Ref. [20| for the main 
results). A convenient formula for our purpose is 
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which describes accurately the first lobe of ji2 for small 
values of b. In order for ji to have an imaginary part the 
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FIG. 1: Imaginary part of the critical exponent \i as function 
of a for b = 0.25. Notice the main lobe centered around a = 1. 
For the first lobe, the difference between the approximation 
(|13p and the numerics is smaller than the thickness of the line. 
The second lobe, shown in the inset, is so much smaller than 
the first one that they can hardly be seen on the same scale. 



argument of arccos must be larger than one (in absolute 
value). Naturally, identifying the most unstable mode 
amounts to finding the extremum value of 



cos 



(ir^a) + 



irb 2 sin (^y^) 
16y/a(a- 1) ' 



To leading order in b this corresponds to a = 1. 

Finally, numerically solving Eq. ([T2]h one readily ob- 
tains 27r/fc, which represents the spacing of adjacent max- 
ima, herein called <S, as a function of the driving fre- 
quency uj of the transverse confinement uj r . Neglecting 
the k 4 term in Eq. ([12]), which is numerically small under 



typical experimental conditions, one has that 
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(2a s p)- 1 / 2 (l + 2a s p) 3 / 4 
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x(4 + 6a s p)- 1 / 2 , 



(14) 



where p = N/2L is the density of the condensate. 

It is important to note that the above formula for k 
[Eq. ([HI) ] has been obtained by assuming a homogeneous 
condensate with constant density p = po- This approxi- 
mation yields a spacing between adjacent maxima of 



a _ 27r 



(15) 



where k is computed using Eq. (p~4|) with p = p$. How- 
ever, the density of the condensate in the considered sys- 
tem is not homogeneous. To account for this inhomo- 
geneity, as a first-order approximation, one uses the fact 
that, typically for the cases under consideration, the den- 
sity of the condensate varies on a space scale much larger 
than that of the observed patterns and thus we can ex- 
tend Eq. ([T4|) with the density being space dependent: 
p = p(x). In fact, the density of the condensate can be 
approximated in this slowly- varying limit by the so-called 
Thomas-Fermi (TF) approximation: 



p{x) 



4L 3 



(16) 



when — L < x < L and p(x) = otherwise. Therefore, it 
is possible to approximate the average spacing by taking 
a spatially averaged k for k in Eq. (fT4|) using p(x) given 
by the TF approximation in Eq. ([T6]) : 
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where 



k = — — j k(x)dx. 
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(17) 
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Alternatively, one could average the spacing directly by 
using the expression 



s 2 = ± 

2 2L 



2tt 



dx, 



(19) 



where, again, k(x) is given by Eq. (Q3j) with the Thomas- 
Fermi density p{x) of Eq. ([16]). 

The above expressions for the spacing S provide the 
analytical prediction of the present study that can be 
readily compared quantitatively with the experimental 
results of Ref. [l2[. This comparison can be seen in 
Fig. [2j where the theoretical predictions for So, Si and 
S2 defined above are adapted to the 87 Rb experiments 
of Ref. [l2[, with ccV 5 o/(27r) = 160.5 Hz, the condensate 
length 2L = 180 pm and TV = 5 x 10 5 atoms [21]. We 
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FIG. 2: (Color online) Average spacing of adjacent maxima 
of the longitudinal patterns as a function of the transverse 
driving frequency. Blue diamond dots depict the experimen- 
tal data of Ref. [12[ (the experimental error bar depicted here 
only takes into account the error bar in the pixel size of the 
recorded experimental image). Dashed black line, solid red 
line and blue dashed-dotted line correspond, respectively, to 
spacings computed using So, Si and S2. Green empty circles 
correspond to spacings extracted from the ID NPSE numer- 
ics by the averaged spacing method. Pink squares correspond 
to the spacings extracted from the full 3D numerics using 
the FFT method and its associated error bars (see text for 
details). Notice that for o;/27r close to 160.5 Hz the theoret- 
ical prediction is far from the experimental data. The radial 
breathing mode excited at these frequencies cannot be cap- 
tured by the NPSE but is well captured by the 3D numerics. 



observe a very good qualitative and a good quantita- 
tive agreement between the theoretical predictions and 
the experimental result, solidifying our expectation that 
Eq. ([T4|) captures accurately and in a fully analytical 
way the observed phenomenology of the experiment of 
Ref. [12]. Is it worth mentioning that Si and S2 are 
closer to the experimental data since they account for the 
inhomogeneity of the cloud. Notice as well that our an- 
alytical prediction shows deviations from the experimen- 
tal data at low frequencies, where the spatial periods of 
the Faraday waves are comparable with the length of the 
condensate. At larger frequencies we have good agree- 
ment between the theoretical curve and the experimental 
data, as the periods of the Faraday waves are substan- 
tially smaller than the spatial extent of the cloud. The 
main sources of slight disparity between the theoretical 
predictions and the experimental results can be traced in 
the transverse Gaussian (as opposed to Thomas-Fermi) 
profile and the fact that the analysis cannot directly in- 
corporate the weak longitudinal trapping potential (see 
the discussion above). These will be further clarified be- 
low, through the comparison with the direct numerical 
simulation results of both the NPSE as well as the full 
3D GP equation. 
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FIG. 3: Faraday pattern formation for the ID NPSE. The 
top panel depicts the growth rate G (L 2 norm of the devia- 
tion from the initial density) of the pattern. The middle two 
rows display the Faraday patterns at the times indicated (see 
vertical lines in top panel) where the left subpanels show the 
density profile while the right subpanels show the deviation 
from the initial density profile. The bottom panel depicts the 
space-time evolution of the density. This case corresponds to 
the experiment in Ref. [l2[, namely, a cloud of N = 5 x 10 5 
87 Rb atoms, trapped by {clV.o, co z }/(2tt) — {160.5, 7} Hz with 
a 20% modulation of the radial confinement at a frequency 
u/(2tt) = 321 Hz. 



III. NUMERICAL RESULTS 

A. One-Dimensional Numerics on the NPSE 

Having completed the linear stability analysis, let us 
now turn to full numerical simulation to investigate the 
instability onset and the emergence of the relevant Fara- 
day patterns. We have simulated the NPSE ([3]) for the 
experimental conditions described in Ref. [l2[. Specifi- 
cally, in Fig. [3l we show the formation of the Faraday 
pattern for a condensed cloud of TV = 5 x 10 5 87 Rb 
atoms contained in a magnetic trap with frequencies 
{^r,o/(27r), (jJ z /(2tt)} = {160.5,7} Hz where the radial 
trap frequency has a modulation of 20% (e = 0.2, which 
is within the typical range of experimentally used modu- 
lations) and a frequency uj/(2tt) = 321 Hz corresponding 
to the resonant oscillation fre quen cy of the radial breath- 
ing mode (i.e., u ~ 2u; r ^) [22ll23j]. As it can be observed 
from the figure, the Faraday pattern grows exponentially 
until it is clearly visible in the density space-time evolu- 
tion (bottom panel) after about 125 ms. It is reassuring 
that the NPSE is successful in capturing the Faraday 
pattern with the same wavelength of 10-11 /im as the 
experiment of Ref. [l2[. On the other hand, we have 



FIG. 4: Same as in Fig. [3] for a (stronger) 40% modula- 
tion of the radial confinement at a (slower) driving frequency 
cj/(2tt) = 150 Hz. 



found that the NPSE cannot capture the right time for 
the development of the instability. The NPSE results 
take about 125 ms (i.e., about 40 periods of the modu- 
lation drive) for the Farad ay p attern to be visible while, 
in the experiments of Ref. [12], some of the patterns are 
clearly visible after some 10 periods of the drive. 

A possible explanation for the discrepancy of the Fara- 
day pattern growth between NPSE and the experiments 
may lie in the size of the initial perturbation, that will 
eventually seed the Faraday pattern. We have used var- 
ious amplitudes for the initial perturbation after we ob- 
tained the steady state solution to Eq. (j4j) by imaginary 
time relaxation. In the results presented in this work we 
used an initial perturbation with an amplitude randomly 
chosen in an interval 0.001 times the local density. We 
also tried larger perturbations, up to ten times larger, 
and the effect is to accelerate the appearance of the pat- 
terns (results not shown here), however we were unable 
to see distinguishable patterns earlier than 30-35 driving 
periods for the above setting. Another effect that needs 
to be taken into account is the amplitude of the modula- 
tion drive. While in the experiments of Ref. [l2[ Faraday 
patterns, for the resonant frequency cj/(2tt) = 321 Hz, 
quickly formed for even small drive amplitudes (less than 
4%, i.e., e < 0.04), our numerical results using the NPSE 
always needed a much larger drive amplitude (e = 0.2 
in the results of Fig. [3]). Therefore, it is clear that the 
NPSE, although clearly able to capture the nature (wave- 
length) of the Faraday pattern, it is unable to predict the 
growth rate of the instability. The reason for this short- 
coming stems from the fact that the transverse compo- 
nent of the wavefunction in Eq. (j3]) is considered to be 
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FIG. 5: (Color online) Method of averaging spacings in the 
inhomogeneous cloud for the NPSE model. The left and right 
columns of plots correspond, respectively, to co/(2ty) = 200 Hz 
and cj/(27t) = 320 Hz. The top row of panels shows the 
density deviation Af from the initial density profile. The 
middle row of panels depicts with dots the measured spacing 
between maxima of Af inside the cloud. The averaged spac- 
ing is shown with the red horizontal line. The bottom row of 
panels shows the FFT of the density |/| 2 . The vertical lines 
depict the estimated window of the frequencies contained in 
the Faraday pattern due to the non-homogeneity in the spac- 
ings (see middle row). 



at its ground state at all times. This is quite a strong 
assumption considering that the cloud presents impact 
oscillator-type dynamics for its radial width [l2[ (cf. top 
panel of Fig. [6] below). This will be verified when we re- 
lax the radial wavefunction profile in the 3D simulations 
shown below. 

The results presented in Fig. [3] correspond to the most 
pattern- forming sensitive case since the drive of the radial 
frequency is tuned to resonate with the natural breath- 
ing frequency of the radial mode. For other driving fre- 
quencies, the growth of the Faraday pattern is less pro- 
nounced. This is demonstrated in Fig. [4j where we use 
the same parameter values as in Fig. [3] but changed the 
driving frequency to u/ (2tt) = 150 Hz and we doubled its 
amplitude (e = 0.4). For this out-of-resonance frequency, 
the Faraday pattern takes longer to form and even a drive 
with twice the amplitude takes longer to seed the pattern 
(the pattern is not visible until approximately 140ms). 



Nonetheless, it is interesting that this out-of-resonance 
case only takes about 20 periods to manifest itself. 

To summarize the results of the ID NPSE simulations 
and to compare them with our analytical prediction for 
the spacing of the ensuing Faraday patterns, we proceed 
to measure the averaged spacing in the simulations. The 
method relies on computing the spacing between max- 
ima on the density deviation Af from the initial profile 
(see top row in Fig. [5]). A couple of examples of the 
dependence of the spacing inside the cloud are depicted 
in Fig. [5] together with their average, from now on de- 
noted as S 7v . It is clear from these examples that at the 
center of the cloud, where the density is larger, the spac- 
ing is larger than at the periphery of the cloud where 
the density is lower. The averaged spacing Sn was com- 
puted as a function of the driving frequency and it is de- 
picted by the green empty circles in Fig. [2l As it is clear 
from Fig. [2j the spacing Si [computed using the analyt- 
ical expression for the spacing given in Eq. ([T4]) with a 
spatially averaged density on the TF approximation] and 
the averaged spacing §n from the NPSE dynamics are in 
good agreement (and the relevant approximation of using 
Eq. (p~4|) together with Eq. ([T6]) to represent the effects 
of the longitudinal potential is a fairly accurate one). It 
is important to mention that the non-homogeneity of the 
spacings induces an inherent uncertainty in the quantifi- 
cation of the associated spacing for a particular Faraday 
pattern. The associated window of spacings contained 
in the Faraday pattern can also be seen from the fast 
Fourier transform (FFT) spectrum of the density. In the 
bottom panels of Fig. [5] we depict the FFT spectrum of 
the density for a couple of cases with their respective 
frequency windows (see below for further elaboration on 
this effect). 



B. Three-Dimensional Numerics 

In order to more accurately model the Faraday pattern 
formation, we used direct numerical simulations of the 
Gross-Pitaevskii equation (pQ). The numerics consist of 
integrating Eq. (pQ) in cylindrical coordinates. The choice 
of cylindrical symmetry instead of full 3D numerics is 
justified by the fact that the radial direction does not de- 
velop azimuthal instabilities as it can be observed in the 
experiments of Ref. [12j] (and it was also verified by addi- 
tional tests runs of the full 3D equation on a coarser grid). 
The main challenge in numerically integrating Eq. (pQ) 
stems from the impact oscillator-type dynamics of the 
radial profile that is driven at resonance. These oscil- 
lations produce two numerically challenging effects: (a) 
they bring most of the atoms close to r = requiring an 
extremely fine grid, and (b) as the cloud accelerates dur- 
ing the impact oscillations, the wavefunction oscillates, in 
space, very rapidly (although the density does not) again 
requiring a very fine grid. Therefore, although (a) could 
be circumvented by a grid refinement around r = 0, chal- 
lenge (b) requires a fine r-grid where the cloud is traveling 
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FIG. 6: Faraday pattern from the (r, z) 3D simulations. 
This case corresponds to the experiment in Ref. [l2j], 
namely, a cloud of N — 5 x 10 5 87 Rb atoms, trapped by 
{oo r ,o,cOz}/(27v) = {160.5, 7} Hz with a 20% modulation of the 
radial confinement at a driving frequency u/(2tt) = 321 Hz. 
The top panel depicts the transverse radius of the cloud dis- 
playing impact-oscillator behavior (thin vertical lines depict 
the times of the snapshots shown in Fig. [TJ. The second 
panel depicts the growth of the Faraday pattern while the 
bottom two rows depict the r-integrated density profiles (left 
subpanels) and their deviation from the initial profile (right 
subpanels) . 




fastest and this happens on a large portion of the domain. 
Thus a fine grid needs to be implemented throughout the 
(radial) r-direction, and as a (numerical scheme stability) 
consequence a very small time step is also required. For 
the (longitudinal) ^-direction, it suffices to have enough 
points to accurately capture the Faraday pattern whose 
wavelength is quite manageable. In the 3D simulations 
shown in Figs. [6] and [71 we were able to accurately inte- 
grate Eq. (pQ) for about 7 cycles of the drive with a grid 
of 2001x401 points in the (r, z)-plane with a finite dif- 
ference scheme in space with 4-5 th Runge Kutta in time 
with a time step of 0.00025 (in adimensionalized units). 

The initial condition used in the simulations (i.e., the 
ground state of the condensate) was obtained, as for the 
ID case, by imaginary time relaxation. 

Figures [6] and [71 depict the Faraday pattern arising 
from the (r, z) Gross-Pitaevskii simulation for a cloud of 
TV = 5 x 10 5 87 Rb atoms driven at resonance (uo/(2tt) = 
321 Hz) by a modulation amplitude of 20%, correspond- 
ing to the experiments of Ref. [l2[. Depicted in Fig. [71 
is the ^/-integrated (top) view which is what is measured 
in the experiments. Clearly observable are the well sep- 
arated fringes of the forming Faraday pattern with an 
approximate spacing of about 8.5 fim. Furthermore, in 
contrast with the NPSE simulations, the Faraday insta- 
bility develops more rapidly (see second panel in Fig. [6j) 
and it is clearly observable after only 6-7 periods of the 




z (urn) 25 



FIG. 7: Development of the Faraday pattern in our 3D numer- 
ical simulations. Shown are the snapshots of the ^/-integrated 
profile density (i.e., the observable in the experiments) for the 
same data as in Fig. [6] at the times indicated. 



drive. As mentioned above, the instability sets in much 
more rapidly in the full 3D system than in the ID NPSE 
reduction, as expected based on the previous discussion. 
It is worth mentioning that in the 3D simulations we did 
not introduce a perturbation to the initial condition to 
seed the Faraday patterns since the inherent numerical 
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FIG. 8: (Color online) FFT analysis for the data from the 
experiments of Ref. [12[. Left and right column correspond, 
respectively, to a driving frequency of uj/(2tt) = 140 Hz and 
uj/(2iz) = 300 Hz. The top row depicts the images from the 
experiment. The second row depicts the corresponding FFTs 
and the third row the FFTs integrated over the y direction. 



noise was capable of starting the pattern. 

To validate our full 3D numerics we have computed 
the average spacing (using the FFT method, see below) 
for three different cases of the driving frequency. The 
results are depicted by the pink squares in Fig. El As it is 
evidenced from the figure, the 3D numerical simulations 
accurately reproduce the spacing from the experiments 
in Ref. [12| including the case when the external driving 
frequency is half of the resonant frequency (see left-most 
pink square point in Fig. [2]). 

In Ref. [IH, the sudden drop in the spacing around 
u/(2tt) = 160 Hz, i.e., half of the resonant frequency, 
was attributed to excitation of the radial breathing mode. 
However, our numerics suggest that this is not the case 
and that this is due to the fact that we are driving a 
sub-harmonic that excites the resonance. 

It is important to note that the data for the experiment 
in Ref. [12[ has a significant variability in the spacing 
values. This natural experimental variability has many 
potential sources. Here, we would like to focus on the er- 
ror generated by the width of the frequency peak used to 
measure the spacing through the fast Fourier transform 
(FFT). In the experiments of Ref. [12[ (see blue diamond 
dots in Fig. [2j), the spacing of the Faraday pattern was 
measured by computing the FFT of the integrated den- 
sity image and extracting the spatial frequency of the 
dominant peak. The error bars depicted in Fig. [2] for the 
experimental data only take into account the error bar in 
the pixel size of the experimental snapshots (and not the 
variability due to the width of the FFT peak, see below). 

In Fig. [8] we present a couple of examples (for uo/ (2tt) = 
140 Hz and uo/ (2tt) = 300 Hz) of the Faraday patterns ob- 
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FIG. 9: (Color online) FFT analysis from the 3D numerics 
for uo/ (2tt) = 321 Hz. The top two panels depict, respectively, 
the Faraday pattern and its FFT. The third panel depicts the 
FFT integrated over the y direction. The bottom panel de- 
picts the same as the third panel but the FFT is computed 
after adding a 5% noise to the density to simulate the experi- 
mental noise in the picture. The vertical lines in the last panel 
represent our estimated window of possible frequencies for the 
Faraday pattern (see error bars in the pink square points in 
Fig.E}. 



served in the experiments of Ref. [12[. Also in the figure, 
we depict the FFT of the data (second row) as well as 
its ^/-integrated counterpart (bottom row). As it is clear 
from the figure, there is a dominant spatial frequency 
that can be measured in order to extract the associated 
spacing of the Faraday pattern. Nonetheless, it is impor- 
tant to note that the dominant peak in the FFT spectrum 
has a width that indicates an interval of spacings that 
make up the original Faraday pattern. The width of this 
peak gives an indication of spatial variability of the pat- 
tern spacing: the wider the spectral peak the more spatial 
variability there exists in the pattern spacing. We have 
observed the same phenomenology when using our 3D 
numerical data. A typical example of the FFT analysis 
of our 3D numerical data is presented in Fig. We have 
checked that our numerical data is able to reproduce the 
behavior of the FFT analysis of the experimental data. 
For a better comparison with experiments we added a 
5% random noise to the numerically computed density 
so as to emulate the noise in the experiment (see bottom 
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panel of Fig. [9]). In order to extract the average spacing, 
S3D , from the 3D numerics we computed the position of 
the dominant peak in the FFT (see pink square points 
in Fig. [2]). As it is clear from the integrated FFT curves 
(see bottom two panels in Fig. [9]), the Faraday pattern 
contains a dominant peak with a respective finite width. 
The width of the peak indicates the presence of a range 
of spatial frequencies (instead of a single one) and thus 
we can associate an error bar to the average spacing by 
taking the width of the dominant peak as the variability 
in spatial frequencies. This variability has been incor- 
porated in the average spacing of 3D numerical data in 
Fig. [2] by means of the vertical error bars for the pink 
square points. In the same spirit, the actual error bar 
in the experiment should be slightly larger to also ac- 
count for variability of the spatial frequency due to the 
width of the computed FFT peaks from the experimental 
data. Within this variability that, based on Figs. [8] and [9l 
clearly exists both in the experimental and the numerical 
(3D) data, we can conclude that our theoretical results 
are in very good agreement with both their experimental 
and their numerical counterparts. 

IV. CONCLUSIONS 

In this communication, we have revisited the topic of 
Faraday waves and corresponding resonances in Bose- 
Einstein condensates. In particular, we have focused 
on the experimentally relevant case where the transverse 
confinement is periodically modulated in time. We have 
used the non-polynomial Schrodinger equation as a tool 
that permits to present in an explicit form this transverse 
modulation in an effective longitudinal equation describ- 
ing the dynamics of the condensate wavefunction. Then, 
a subsequent modulational stability analysis permitted 
us to examine the stability of spatially uniform states in 
this transversely driven, yet effectively one-dimensional 
setting. This analysis, leading to a Mathieu equation, 
combined with the well-established theory of the latter 
equation allowed us to identify the dominant mode of the 
instability. This, in turn, led us to extract an explicit an- 
alytical formula that allows for this mode's wavenumber 
(and hence its wavelength which is directly associated 
with the pattern periodicity and hence is an experimen- 
tally observable quantity) as a function of the driving fre- 
quency of the transverse confining potential. Direct com- 
parison of the fully analytical result with the experimen- 



tal observations confirmed the accuracy of our approach. 
These analytical and experimental results were also cor- 
roborated by numerical computations both within the 
framework of the one-dimensional NPSE equation, as 
well as for the case of the fully 3D Gross-Pitaevskii equa- 
tion. The similarities of the two regarding the instability 
length scale and the differences of the two in connec- 
tion to the instability growth rate have been accordingly 
highlighted. These computations have allowed us to val- 
idate the quality of our theoretical approximations and 
give a detailed comparison between theory, numerics and 
experiment. 

There are numerous interesting possibilities that this 
experiment presents for future studies. One of them 
is to consider the predominantly two-dimensional case 
of pancake-shaped condensates, where, depending on 
the driving frequency, square or rhombic patterns may 
form. In that setting too, the effective wave equations of 
Ref. [ItJ may allow to perform the modulational stability 
analysis and obtain a quantitative handle on the domi- 
nant unstable mode. On the other hand, modulations 
of all three spatial dimensions of the confining potential 
would be of interest in their own right. In the latter 
case, while reductions of the type used herein would not 
be relevant, nevertheless the dynamics may still be ana- 
lytically describable upon appropriate assumptions, such 
as, e.g., the quadratic phase assumption of Ref. [5|, by 
means of coupled, nonlinear ordinary differential equa- 
tions such as Eqs. (12) of Ref. @. Understanding these 
settings in more quantitative detail, both analytically, 
numerically and experimentally, is under current exami- 
nation and will be reported in future publications. 
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